library(dplyr)
#library(psych) #for pairs.panels, but could use other packages, e.g. GGalley
library(lavaan)
library(semPlot)
library(DiagrammeR)
library(tidyr)
library(ggplot2)

Import data

combined=read.csv("data/annual_averages/annual_data_compiled_regions.csv",stringsAsFactors = F)
cnames=read.csv("analysis/column_names_region.csv", stringsAsFactors = F)
dsub=filter(combined, Year>=1975) %>% arrange(Region,Year)
focaldata=dsub[,cnames$Datacolumn]
fvars=cnames$Shortname
colnames(focaldata)=fvars
regions=unique(focaldata$region)
regionorder=c("West","North","South")
years=1975:2021

focaldata$tzoop=focaldata$hcope+focaldata$clad+focaldata$mysid+focaldata$pcope
fvars=c(fvars,"tzoop")
cnames=rbind(cnames,data.frame(Longname=NA,Shortname="tzoop",
                               Diagramname="Total Zooplankton\nBiomass",
                               Datacolumn=NA,Log="yes"))

#focal variables
varnames=c("temp", "flow","nitrate","ammonia","dophos","chla","hcope","clad","amphi","pcope","mysid","potam","corbic","sside","estfish_bsot","estfish_bsmt","tzoop")

source("analysis/myLavaanPlot.r")

Data prep

Log transform, scale

#log transform
logvars=fvars[cnames$Log=="yes"]
logtrans=function(x) {
  x2=x[which(!is.na(x))]
  if(any(x2==0)) {log(x+min(x2[which(x2>0)],na.rm=T))}
  else {log(x)}
}
focaldatalog = focaldata %>% 
  mutate_at(logvars,logtrans)

#scale data
fdr0=focaldatalog
tvars=fvars[-(1:2)]

fdr=fdr0 %>% group_by(region) %>% 
  #lag
  mutate_at(tvars,list("1"=lag)) %>% 
  #scale
  mutate_at(-(1:2),scale) %>% 
  ungroup() %>% 
  as.data.frame() 

#detrended data
fdr_dtr=fdr0 %>% group_by(region) %>% 
  #detrend
  mutate_at(tvars,function(x) { 
    x<<-x
    if(!all(is.na(x))) {
      if((length(which(x==0))/length(x))<0.5) {
        x2<<-x
        x2[x2==0]=NA
        res<<-residuals(lm(x2~years))
        out=x
        out[which(!is.na(x2))]=res
        return(out)
      } else {return(x)}
    } else {return(x)}
  }) %>%
  #lag
  mutate_at(tvars,list("1"=lag)) %>% 
  #scale
  mutate_at(-(1:2),scale) %>% 
  ungroup() %>% 
  as.data.frame() 

Time series plots

Cross-correlation matrices

(only sig correlations shown… no correction for multiple comparisons)

Fish indices are not correlated in S and N!

SEM model

With and without detrending.

#west has no ssides, corbic

# modwest='zoop=~hcope+mysid
#         fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+potam+flow
#         chla~potam+flow
#         fish~zoop+flow
# '

modwest='chla~potam+flow
        tzoop~chla+potam+flow
        estfish_bsmt~tzoop+flow
        estfish_bsot~tzoop+flow
'

modfitwest=sem(modwest, data=filter(fdr,region=="West"))
modfitwest_dtr=sem(modwest, data=filter(fdr_dtr,region=="West"))
summary(modfitwest, standardized=T, rsq=T)
## lavaan 0.6-4 ended normally after 27 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         14
## 
##                                                   Used       Total
##   Number of observations                            40          47
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                       4.813
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.307
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     potam            -0.325    0.135   -2.407    0.016   -0.325   -0.380
##     flow              0.098    0.131    0.748    0.454    0.098    0.118
##   tzoop ~                                                               
##     chla              0.831    0.096    8.676    0.000    0.831    0.792
##     potam            -0.154    0.088   -1.756    0.079   -0.154   -0.171
##     flow             -0.069    0.080   -0.867    0.386   -0.069   -0.080
##   estfish_bsmt ~                                                        
##     tzoop             0.626    0.115    5.457    0.000    0.626    0.568
##     flow              0.386    0.100    3.874    0.000    0.386    0.403
##   estfish_bsot ~                                                        
##     tzoop             0.598    0.119    5.030    0.000    0.598    0.542
##     flow              0.389    0.103    3.771    0.000    0.389    0.407
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_bsmt ~~                                                       
##    .estfish_bsot      0.054    0.066    0.811    0.418    0.054    0.129
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .chla              0.584    0.131    4.472    0.000    0.584    0.802
##    .tzoop             0.215    0.048    4.472    0.000    0.215    0.267
##    .estfish_bsmt      0.402    0.090    4.472    0.000    0.402    0.412
##    .estfish_bsot      0.430    0.096    4.472    0.000    0.430    0.442
## 
## R-Square:
##                    Estimate
##     chla              0.198
##     tzoop             0.733
##     estfish_bsmt      0.588
##     estfish_bsot      0.558
summary(modfitwest_dtr, standardized=T, rsq=T)
## lavaan 0.6-4 ended normally after 24 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         14
## 
##                                                   Used       Total
##   Number of observations                            40          47
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                       6.849
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.144
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     potam             0.032    0.162    0.196    0.844    0.032    0.034
##     flow              0.215    0.160    1.341    0.180    0.215    0.229
##   tzoop ~                                                               
##     chla              0.837    0.105    8.009    0.000    0.837    0.795
##     potam             0.026    0.107    0.247    0.805    0.026    0.026
##     flow             -0.020    0.108   -0.183    0.854   -0.020   -0.020
##   estfish_bsmt ~                                                        
##     tzoop             0.417    0.115    3.618    0.000    0.417    0.428
##     flow              0.446    0.114    3.915    0.000    0.446    0.463
##   estfish_bsot ~                                                        
##     tzoop             0.379    0.121    3.139    0.002    0.379    0.383
##     flow              0.455    0.119    3.824    0.000    0.455    0.467
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_bsmt ~~                                                       
##    .estfish_bsot      0.046    0.089    0.519    0.604    0.046    0.082
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .chla              0.879    0.197    4.472    0.000    0.879    0.953
##    .tzoop             0.384    0.086    4.472    0.000    0.384    0.376
##    .estfish_bsmt      0.534    0.119    4.472    0.000    0.534    0.548
##    .estfish_bsot      0.584    0.131    4.472    0.000    0.584    0.585
## 
## R-Square:
##                    Estimate
##     chla              0.047
##     tzoop             0.624
##     estfish_bsmt      0.452
##     estfish_bsot      0.415
# par(mfrow=c(1,2))
# semPaths(modfitwest, "std", edge.label.cex = 1, residuals = F)
# semPaths(modfitwest, "par", edge.label.cex = 1, residuals = F)

labelswest <- createLabels(modfitwest, cnames)

# residuals(modfitwest)
# modificationindices(modfitwest)
# modnorth='zoop=~hcope+mysid
#         #fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+potam+flow
#         chla~potam+flow
#         estfish_bsmt~zoop+flow
#         estfish_bsot~zoop+flow
# '
# modnorth='zoop=~clad
#         zoop~chla+corbic+potam+flow
#         chla~corbic+potam+flow
#         estfish_bsmt~zoop+flow+sside+chla
#         estfish_bsot~zoop+flow+sside+chla
# '

modnorth='chla~corbic+potam+flow
        tzoop~chla+corbic+potam+flow
        estfish_bsmt~tzoop+flow+chla
        estfish_bsot~tzoop+flow
'

modfitnorth=sem(modnorth, data=filter(fdr,region=="North"))
modfitnorth_dtr=sem(modnorth, data=filter(fdr_dtr,region=="North"))
summary(modfitnorth, standardized=T, rsq=T)
## lavaan 0.6-4 ended normally after 22 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         17
## 
##                                                   Used       Total
##   Number of observations                            40          47
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                       1.857
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.868
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.358    0.132    2.718    0.007    0.358    0.422
##     potam            -0.121    0.143   -0.844    0.399   -0.121   -0.140
##     flow              0.022    0.138    0.159    0.874    0.022    0.026
##   tzoop ~                                                               
##     chla              0.681    0.112    6.072    0.000    0.681    0.614
##     corbic            0.414    0.102    4.069    0.000    0.414    0.440
##     potam             0.075    0.103    0.735    0.463    0.075    0.078
##     flow             -0.351    0.098   -3.586    0.000   -0.351   -0.370
##   estfish_bsmt ~                                                        
##     tzoop             0.033    0.188    0.175    0.861    0.033    0.032
##     flow              0.040    0.130    0.304    0.761    0.040    0.041
##     chla              0.694    0.215    3.227    0.001    0.694    0.614
##   estfish_bsot ~                                                        
##     tzoop             0.468    0.135    3.476    0.001    0.468    0.470
##     flow              0.251    0.128    1.968    0.049    0.251    0.266
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_bsmt ~~                                                       
##    .estfish_bsot      0.025    0.098    0.256    0.798    0.025    0.040
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .chla              0.562    0.126    4.472    0.000    0.562    0.737
##    .tzoop             0.283    0.063    4.472    0.000    0.283    0.301
##    .estfish_bsmt      0.565    0.126    4.472    0.000    0.565    0.579
##    .estfish_bsot      0.677    0.151    4.472    0.000    0.677    0.727
## 
## R-Square:
##                    Estimate
##     chla              0.263
##     tzoop             0.699
##     estfish_bsmt      0.421
##     estfish_bsot      0.273
summary(modfitnorth_dtr, standardized=T, rsq=T)
## lavaan 0.6-4 ended normally after 23 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         17
## 
##                                                   Used       Total
##   Number of observations                            40          47
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                       2.521
##   Degrees of freedom                                 5
##   P-value (Chi-square)                           0.773
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.315    0.162    1.940    0.052    0.315    0.311
##     potam             0.066    0.153    0.431    0.667    0.066    0.070
##     flow              0.105    0.169    0.619    0.536    0.105    0.108
##   tzoop ~                                                               
##     chla              0.710    0.097    7.337    0.000    0.710    0.679
##     corbic            0.317    0.104    3.049    0.002    0.317    0.299
##     potam             0.086    0.094    0.920    0.358    0.086    0.087
##     flow             -0.460    0.104   -4.414    0.000   -0.460   -0.454
##   estfish_bsmt ~                                                        
##     tzoop             0.356    0.181    1.972    0.049    0.356    0.373
##     flow              0.204    0.135    1.507    0.132    0.204    0.211
##     chla              0.333    0.187    1.784    0.074    0.333    0.334
##   estfish_bsot ~                                                        
##     tzoop             0.228    0.149    1.528    0.126    0.228    0.238
##     flow              0.225    0.151    1.492    0.136    0.225    0.232
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_bsmt ~~                                                       
##    .estfish_bsot      0.085    0.109    0.776    0.438    0.085    0.124
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .chla              0.846    0.189    4.472    0.000    0.846    0.871
##    .tzoop             0.317    0.071    4.472    0.000    0.317    0.298
##    .estfish_bsmt      0.526    0.118    4.472    0.000    0.526    0.542
##    .estfish_bsot      0.892    0.199    4.472    0.000    0.892    0.916
## 
## R-Square:
##                    Estimate
##     chla              0.129
##     tzoop             0.702
##     estfish_bsmt      0.458
##     estfish_bsot      0.084
# par(mfrow=c(1,2))
# semPaths(modfitnorth, "std", edge.label.cex = 1, residuals = F)
# semPaths(modfitnorth, "par", edge.label.cex = 1, residuals = F)

labelsnorth <- createLabels(modfitnorth, cnames)

# residuals(modfitnorth)
# modificationindices(modfitnorth)
#no potam
# modsouth='zoop=~hcope+mysid
#         #fish=~estfish_bsmt+estfish_bsot
#         zoop~chla+corbic+flow
#         chla~corbic+flow
#         estfish_bsmt~zoop+flow
#         estfish_bsot~zoop+flow
# '

modsouth='chla~corbic+flow
        tzoop~chla+corbic+flow
        estfish_bsmt~tzoop+flow+corbic
        estfish_bsot~tzoop+flow+corbic
'

modfitsouth=sem(modsouth, data=filter(fdr,region=="South"))
modfitsouth_dtr=sem(modsouth, data=filter(fdr_dtr,region=="South"))
summary(modfitsouth, standardized=T, rsq=T)
## lavaan 0.6-4 ended normally after 20 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         16
## 
##                                                   Used       Total
##   Number of observations                            40          47
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                       1.175
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.556
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.446    0.154    2.902    0.004    0.446    0.436
##     flow             -0.014    0.147   -0.093    0.926   -0.014   -0.014
##   tzoop ~                                                               
##     chla              0.677    0.095    7.107    0.000    0.677    0.705
##     corbic            0.242    0.102    2.371    0.018    0.242    0.246
##     flow             -0.138    0.088   -1.559    0.119   -0.138   -0.147
##   estfish_bsmt ~                                                        
##     tzoop            -0.008    0.164   -0.052    0.959   -0.008   -0.008
##     flow             -0.020    0.140   -0.142    0.887   -0.020   -0.021
##     corbic            0.491    0.170    2.896    0.004    0.491    0.496
##   estfish_bsot ~                                                        
##     tzoop             0.184    0.150    1.231    0.218    0.184    0.191
##     flow             -0.175    0.128   -1.366    0.172   -0.175   -0.193
##     corbic            0.423    0.155    2.727    0.006    0.423    0.446
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_bsmt ~~                                                       
##    .estfish_bsot      0.004    0.108    0.035    0.972    0.004    0.006
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .chla              0.846    0.189    4.472    0.000    0.846    0.813
##    .tzoop             0.307    0.069    4.472    0.000    0.307    0.320
##    .estfish_bsmt      0.745    0.167    4.472    0.000    0.745    0.765
##    .estfish_bsot      0.623    0.139    4.472    0.000    0.623    0.697
## 
## R-Square:
##                    Estimate
##     chla              0.187
##     tzoop             0.680
##     estfish_bsmt      0.235
##     estfish_bsot      0.303
summary(modfitsouth_dtr, standardized=T, rsq=T)
## lavaan 0.6-4 ended normally after 15 iterations
## 
##   Optimization method                           NLMINB
##   Number of free parameters                         16
## 
##                                                   Used       Total
##   Number of observations                            40          47
## 
##   Estimator                                         ML
##   Model Fit Test Statistic                       0.825
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.662
## 
## Parameter Estimates:
## 
##   Information                                 Expected
##   Information saturated (h1) model          Structured
##   Standard Errors                             Standard
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   chla ~                                                                
##     corbic            0.160    0.164    0.973    0.330    0.160    0.157
##     flow             -0.032    0.158   -0.206    0.837   -0.032   -0.033
##   tzoop ~                                                               
##     chla              0.639    0.112    5.708    0.000    0.639    0.642
##     corbic            0.200    0.117    1.707    0.088    0.200    0.198
##     flow             -0.188    0.112   -1.685    0.092   -0.188   -0.193
##   estfish_bsmt ~                                                        
##     tzoop            -0.223    0.153   -1.460    0.144   -0.223   -0.229
##     flow             -0.084    0.149   -0.566    0.571   -0.084   -0.089
##     corbic            0.343    0.158    2.168    0.030    0.343    0.347
##   estfish_bsot ~                                                        
##     tzoop             0.082    0.146    0.562    0.574    0.082    0.085
##     flow             -0.220    0.142   -1.549    0.121   -0.220   -0.236
##     corbic            0.349    0.151    2.316    0.021    0.349    0.360
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .estfish_bsmt ~~                                                       
##    .estfish_bsot     -0.056    0.130   -0.429    0.668   -0.056   -0.068
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .chla              1.010    0.226    4.472    0.000    1.010    0.977
##    .tzoop             0.506    0.113    4.472    0.000    0.506    0.494
##    .estfish_bsmt      0.857    0.192    4.472    0.000    0.857    0.879
##    .estfish_bsot      0.781    0.175    4.472    0.000    0.781    0.829
## 
## R-Square:
##                    Estimate
##     chla              0.023
##     tzoop             0.506
##     estfish_bsmt      0.121
##     estfish_bsot      0.171
# par(mfrow=c(1,2))
# semPaths(modfitsouth, "std", edge.label.cex = 1, residuals = F)
# semPaths(modfitsouth, "par", edge.label.cex = 1, residuals = F)

labelssouth <- createLabels(modfitsouth, cnames)

# residuals(modfitsouth)
# modificationindices(modfitsouth)

Nice plots

Without covariances

Original units

West

North

South

Detrended

West

North

South

With covariances

Original units

Detrended